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Abs tract 

Stark broadening calculations of isolated neutral atom lines in the 
classical path approximation usually involve an electron (moving along a 
straight line path) interacting directly with the atom. Correlations be- 
tween the electrons are then taken into account by imposing a cutoff in 
the interaction, when the distance from the atom, p, exceeds the Debye 
length, p^. A more consistent procedure for the correlation effects is 
to replace the direct interaction of the electron by a Debye-shielded in- 
teraction. The functions A, B, a and b which are required in the theory 
when the Debye-shielded interaction is used are considered in detail in 
this report. When p/p d < 0.1, it is found that the a and b functions may 
be closely reproduced by using unshielded functions in conjunction with 
an upper cutoff of 0.68 p^. In the appendix is a computer program to 
generate these functions written by U. Palmer of JILA, 
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1. Introduction 

Calculations of the broadening by electrons in the impact approxima- 
tion usually involve the evaluation of an operator (the <j> operator) in 
which the classical path S-matrices for the electron collisions have been 
expanded to second order (see, for example, Refs. 1, 2 and 3). Higher 
order terms in the expansion are then accounted for by a lower cutoff in 
the integral over impact parameters, and the effect of electron correla- 
tions is included by an upper cutoff at an impact parameter at the order 
of the Debye length. The purpose of this report is to give an alterna- 
tive numerical procedure to the use of the upper cutoff (although we will 
show in Sec. 6 that a judicious choice of upper cutoff will, under some 

circumstances, reproduce our results. 

4 5 

More recent theory * has shown that the broadening operator for 
the complete line profile which is correct to second order in the inter- 
action potential V(t) can be expressed (for |a>, |a’> matrix elements), 
using Eq. (49) of Ref. 4, as 

r°° iAaj, t 

<a |^(Aw) |a'> = -i£ <aj dte {V(t) Jk><k|v(0) } |a’> . (1.1) 

k J o AV 

Here intermediate states |k> have been explicitly inserted and Aco^. is the 
frequency difference of the radiation from the, intermediate state. This 

\ g j 

result was originally obtained by Baranger from a wing expansion, rather 
than the more complete general theory. In particular this result has im- 
portant consequences , since it shows [see Eq. (1.3) below] that, in the 
correct second-order treatment, functions like A(z is z 2 ) as introduced in 
Ref. 3 corresponding to off-diagonal elements of the broadening operator 
are needed only for the simple case of z^=z^. 
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Putting the potential V(t) equal to the dipole interaction -J a ?(t) 

where E(t) is the electric field due to all the electrons , we see that 
Eq. (1.1) involves the evaluation of the electric field autocorrelation 

function, namely : 

,00 

g(Aco) = dte lAtot {E(t)l(0)> Av (1.2) 

•*o 


This function has been examined in Refs. 7 and 8. In particular, it is 
possible to evaluate the electric field average as if each electron were an 
independent quasi-particle interacting with the radiating atom through its 
dynamically screened electric field. The full evaluation of the dynamic- 
ally screened potential requires the use of the wave number and frequency- 
dependent dielectric constant e + (k,k»v), however, the dominant contribution 
to g(Aw) is from the region where Aw ^ k«v [and in fact* k*v is set equal 

g 

to Aw for the real part of g(Aw) ]* In general , the full evaluation of 

8 9 

g(Aw) is quite complicated, 5 but, when Aw << w^ (the plasma frequency) 
the dynamic dielectric constant can be approximated by the static dielectric 


constant and the shielded field is the Debye-screened field, 
unshielded fields must be used. 


When Aw >> w 9 
p 


In this report, we consider the evaluation of the second-order terms 
using static Debye-screened fields. This is therefore only strictly cor- 
rect when Aw (the frequency separation to the intermediate state) is less 
than w^$ however, correlation and shielding are most important in this region 
where the lines of the spectrum are overlapping or starting to overlap * In 

addition, when Aw >> w the difference between using static Debye-screened 

P > 

fields and unshielded fields in Eq* (1*2) is negligible, and the value of 
Eq, (1,2) is quite small in any case. Finally, by comparison with the re- 
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sults of Ref. 8 [in particular for the real part of g(Aco)] errors in the 


region Aw “p f rom the use °f Debye-screened fields are not expected to be 
large, provided g(Aw) is not appreciably enhanced in this region due to 
instabilities and other non-thermal effects. 

To relate Eq. (1.2) to the usual A and B functions, it is necessary 
to rewrite it slightly. Since the electrons act as independent quasi-part- 
icles the average in Eq. (1.2) may be written in terms of integrals over 
the frequency of collisions dv (an integral essentially over velocities 
and impact parameters) and the time of closest approach t (see for exam- 
ple section 4B of Ref. 4 and Ref. 10). 

Thus 


e lAt0t {E(t)E(0)> Av dt 


dt 


dv 


J 


<■<» where "6 represents the 

dte lAa)t E (t+t )E (t ) ® . 1 ,8 

s o s o shielded field 
o 




■00 

= 

dv 

dx 

z 



— OO 



fiCO 

— 

dv 

dx.. 

J 

J 

00 1 J 


iAfa>(x--x 9 ) x-«t+t 

dx-e E (x-)E (x ) with 

1 si s 2 x 0 -t 

* z o 


x^ iA 
dx 2 e 


by using the Dirichlet integral formula [see Eq. (57), Ref. 4] 


(1.3) 


2. The Shielded Broadening Functions 

Thus, using static Debye-screened fields instead of pure Coulomb 

fields, which amounts to multiplying the latter by (l+r/p^)exp(-r/P D ) , the 

1-3 

usual second-order time integral is replaced by 
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F(z , z ' ;q) = A + iB = 


1 

2 


r 


X . 


dx dx’ (1+xx ' ) g (x ; q) g (x ' ;q)e 


i(zx-z’x') 


, ( 2 . 1 ) 


where 


g(x;q) = l+q(l+x 2 ) 1/2 " (l+x 2 )" 3/2 e' q(1+X } 


2 , 1/2 


( 2 . 2 ) 


2 1/2 

z = pw/v, z' = pu)’/v, q = p/p d> P is the impact parameter, P D = (kT/4irNe ) ' 

the Debye length, and v the electron velocity. It is easily seen that the 

new A and B functions possess the same symmetry properties as the unshielded 
. 3 

ones, vxz. 


A(z,z’;q) = A(z’,z;q) = A(-z,-z';q) 


(2.3) 


B(z,z’;q) = B(z',z;q) - -B(-z,-z';q) 


(2.4) 


The general expression in Eq. (2.1) should be used in the impact 
1 3 

theory ’ to compute the matrix elements of the broadening operator for 
lines with forbidden components (for an isolated line z' = z) . However , 
as stressed earlier, in the unified theory approach only "diagonal" func- 
tions with z = z ' = pAm ^/ v will occur . Of course, no upper cutoff on im- 
pact parameters is required with the new functions. 

Sticking to real integration variables in Eq. (2.1), simple parity 
considerations show that the A-function is given by 


dx g(x;q)cos(zx) dx* g(x ! ;q) cos (z 'x' ) 
o * o 


+ dx xg(x;q)sin(zx) I dx' x'g(x' ;q)sin(z'x') . (2.5) 
o ■’ o 


r 


The necessary Fourier transforms are readily obtained from 


11 


r te( a Wr in e-« aW > 1/2 c<,3M = y.(,V) 1,! l 


( 2 . 6 ) 
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f fc - a«( I V)- 1 / Z Jl |c t ( ! V) 1 «], (2.7) 


' o 

by differentiation with respect to a. Here Kq and are modified Bessel 
functions of the second kind. Hence: 


A(z,z';q) = zz’KqOOKqCr’) + HR'K^RU^R') 

where we have put for brevity 

, 2 2 . 1/2 , , ,2 2 . 1/2 

R = (z +q ) , R' = (z* +q ) ' 


( 2 . 8 ) 


(2.9) 


For the diagonal function we get 

A(z,z;q) = A(z;q) = z 2 K q 2 (R) + R 2 K 1 2 (R) (2.10) 

Since z s z\ and q are all proportional to p s the integration over the 
impact parameter (now from p^ n to °°) can be expressed in terms of 

*00 

a(z,z ? ;q) - j A(Az,Xz ? ;Aq) , (2,11) 

which represents the natural generalization of the function a(z) of GBKO, 

To save space 9 we give here the results only for the case z = z f : 

a (z ;q) = RK Q (R)K ;L (R) - |q 2 [ K ± 2 (R)-K Q 2 (R)] (2.12) 

It is also convenient to introduce in Eqs. (2.10), (2.12), a new 

2 1/2 

parameter z^ = z/q = ppAto/v, so that R = gz, with g “(l+l/z^ ) . For 

z D 'v Aco/tUp >> 1 one has A(z;z/z D ) & A(z) and 3l(z',z/z^) k a(z). Therefore, 
when Ato » tOp, using either shielded fields or unshielded fields in Eq. (1.2) 
gives the same results. As we have said before, the fact that we get the 
correct answer in the important limit of Ato << tOp and that both shielded and 
unshielded fields give nearly identical results when Ato >> is the main 
justification for the functions presented here. 
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3. Dispersion Relations for Off-Diagonal Broadening Functions 
A convenient way to calculate the B function, once A is known, is pro- 
vided by the dispersion relation, which expresses B as a Hilbert transform 
of A, This method has been extensively used in the past to compute various 

diagonal B functions contributing to the broadening of isolated lines emit- 

l 

ted by neutral atoms and positive ions. 5 In both cases the disper- 

sion relations merely reflect the analyticity of the second-order time in- 
tegral with respect to one of its parameters, considered as a complex vari- 
able a 

The extension to off-diagonal functions must be done with some care. 

In this case it is better to introduce an auxiliary complex variable, while 
keeping all physical parameters real (the dispersion relations given re- 
cently"^ for the neutral functions A(z,pz) and B(z,pz), where p = z'/z is 
a real fixed ratio, are wrong, since for p 4 - 1 the time integral does not 
possess the required analytic properties when z is allowed to take complex 
values). For instance, from the structure of Eq. (2.1): 


F(z,z' ;q) 


dx 


dx’ 


f (x,x ? ;q) 


i(zx-z'x') 

e 


(3.1) 


it is apparent that for any real z, z', and q >_ 0, the function <f(0 = 
F(z+£,z'+£;q) is holomorphic in the upper half-plane lm?>0. Applying the 
usual analysis , based upon Cauchy's theorem and the well-known relation 
l/(C-iO) = P(l/S) + iir<5(£), one readily gets 


<?( 0 ) 



dC 


no 


— oo 


c 


9 


(3.2) 
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whence 


B(z s z’ ;q) = - — P 


dC 


A(z+?,z f +?;q) 


(3.3) 


or alternatively 

fCQ 

B(z,z’ ;q) = - ^ P ^(z+p )' AtS 4 ^ 2-2 ') ,C-^(z-z ? ) ;q] • (3.4) 

J —00 

Notice that Eq. (3.3) can be given also the equivalent form 


B(z,z';q) = 


dC 


tA(z+C,z , +?;q)-A(z-?,z , -C;q)] 


(3.5) 


from which the singularity has been removed. 

Although the above relations allow us in principle to compute B to the 
desired accuracy, in practice the limiting process inherent to all of them 
might become sometimes a serious source of trouble. 


4. Complex Integration 

We shall now apply to the shielded case the powerful contour integra- 
tion method which led to closed-form expressions for both A and B in the 

16 

non-shielding limit q = 0. To this end we assume z >_ z' > 0, and make 
the change of variables x = sinh u, x* = sinh u'. Equation (2.1) may then 
be rewritten as 


F = F 1 (z,z';q) + F 2 (z,z';q) , 


(4.1) 


where 


F = — 


r°° 


j . , 1+q cosh u -q cosh u +iz sinh u , , . , , 

du sinh u — a 2 e (Mu) , (4,2) 

cosh u 


with 



ru 


^(u) = 


du * sinh u 


, l+q cosh u’ -q cosh u f -iz ' sinh 

,2 » 6 

cosh u 


and 


F = — 
S 2 2 


, l+q cosh u -q cosh u + iz sinh u , , N 
du — 3 2 e <f> 2 ^ 

-CO cosh u 


with 


+ 2 (u) “ 


u 


, l+q cosh u' e ~q cosh u' - iz'sinh u' 
cosh^u' 


Let us first evaluate F^. Since 


sinh x cosh % e ~q cosh x 
cosh^x 


d / e ~ q cosh x \ 

dx \ cosh x ) * 


integration by parts gives 


, \ . -q cosh u - iz'sinh u . , , , v 

$^(u) = -(1/cosh u) e - iz' ij^(u) , 


where 


, z « f U , , -q cosh u' - iz' sinh u' 
^(u) = du' e M 


and further 


2F 1 = - 


j , l+q cosh u -2q cosh u + iq sinh u 
du sinh u — ^ ^ e 

cosh u 


fiCO 

. . ... , N -2q cosh u + iri sinh u . 

- iz' (du/cosh u) e M + zz'G^ 

* —CO 


where n = z-z' > 0, and 


G 1 " 


du e 


, , . . , ,u , , -q cosh u' - iz’ 

-q cosh u + xz sinh u du e 


U? ,( 4 . 3 ) 


, ( 4 . 4 ) 


( 4 . 5 ) 


( 4 . 6 ) 


( 4 . 7 ) 


( 4 . 8 ) 


( 4 . 9 ) 


sinh u' 


—oo 


—00 


( 4 . 10 ) 





Figure 1 
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fU 


segment CD vanishes , since 
fore 


du 5 ... is an odd function of u, and there- 


du e 


-R cosh u 


du'e 


, -R* cosh u’ 


rn-a-a 


9 


+ i d0 e 


-R ? cosh (u+i0) 


= 2 K 0 (R)K 0 (R f ) + i Mj_ 


(4,13) 


where 


fiT-a-a 


M i = 


r°° 


de 


du e 


-R cosh u - R’ cosh(u+i0) 


(4.14) 


After shifting the u integration from the real axis onto the line 
Im u = - X( u u-ix) , with x given by tan x = R'sin 0/(R + R'cos0) > 0, 
and changing 0 in it - 20, one gets 

nr/2 

\ - 4 I d0 K Q [ S ( 0 ) ] , (4.15) 

(a+a 5 ) /2 

where 


see) = (R 2 + R ’ 2 


- 2RR’ cos 20)' 


(4.16) 


The second term in Eq* (4*1) can be treated in a similar manner , 
but the calculation is slightly more cumbersome * Starting from 


1+q cosh x ^-q cosh x 
cosh^x 



sinh x -q 

0 * 

cosh x 


cosh xl 


+ q cosh x e 


~q cosh x 


(4.17) 
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one first obtains by partial integration 


2 f 2 = 


f°° j , , 1+q cosh u “2q cosh u + in sinh u 

du exnh u — 3 ^ e 

cosh u 


(4.18) 


r 


, ^ , / . , . i , , v -2q cosh u 4- in sinh u , 

du tanh u (q cosh u 4* iz sinh u) e 4* RR ? G 9 , 


where 


G 2 ~ 


du sinh (u4ia) e 


iR sinh (u4ia) 


ru 


. , . , / ? , ,v -iR ? sinh(u f -ia T ) 

du* sinh(u f -ia' )e , 

(4.19) 


By contour integration this can be reduced further to 


G 2 = 2 K 1 (R)R 1 (R ! ) 4 i v, 


(4.20) 


where 


M 2 - 4 


7t/2 
( a + a ' )/2 


d0 sin 0 K q (S) + 4 


rir/2 

( a + a ' ) /2 


d0 K (S)/S 


-4(R-R') 2 


p/2 


de cos e K 2 (s)/s i 


(4.21) 


( a + a ' )/2 


With. F^ determined from Eqs. (4.9), (4.13) and (4.15), and F 2 from Eqs. 
(4.18), (4.20) and (4.21), we go back to Eq. (4.1) and separate the real and 
imaginary parts of F. This yields for A the closed form expression already 
found in Sec. 2: 

A(z ,z ' ;q) = z z , K 0 (R)K 0 (R') + R R'K^RjK^R') , (4.22) 


and for B the integral representation 
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B(z,z’;q) = -q(z+z ' )K 1 (S Q ) /S Q + 2zz* 


"tt/2 
(a+a' ) /2 


d6 K Q (S) 


+ 2RR’ 


fTt/2 


de sin e K q (S) + 2RR’ 


(a+a’)/2 


- 2RR’(R-R') 2 


n r/2 


de cos e k 2 (s)/s 


tt/2 

(a+’a* ) /2 
2 


d6 ^(S) /S 


(4.23) 


(a+a ’ ) /2 


2 2 ^ 

where Sq = (n +4q ) 2 , and S is given by Eq. (4.16). 

The above form of B is not quite satisfactory for numerical computa- 
tion 9 since it contains terms which diverge individually under certain cir- 
cumstances 5 although the net result will always be finite (for instance, in 
the diagonal case z 1 -z the first term behaves like 1/q when q 0,. and so 
does the fourth)* This difficulty is, however, easily resolved by noticing 

that Sq is nothing but the value of S at the lower limit of integration 
0 = %(a+a f ). Developing the right-hand side of the identity 


-q (z+z ? )K^(Sq) /Sq 



(a+a f ) /2 


d 


sin 20 



(4.24) 


and substituting into Eq a (4.23) then enables us to eliminate the unpleasant 
terms and eventually results in the new integral representation 


B(z,z' ;q) 


f tt/2 


2zz 


d0 K q (S) - 2RR’ 


(a+a')/2 


• tt/2 

d6 cos 26 K„(S) , (4.25) 
J u 

(a+a ’ ) /2 


which is the best we can get. 

In view of the subsequent integration over the impact parameter we de- 
fine, by analogy with Eq. (2.11), the function 
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b (z ,z ' ;q) 


dX 


B(Xz,Xz' ;Xq) 


(4.26) 


which generalizes the function b(z) of GBKO. Substituting from Eq. (4.25) 
and noticing that a and a' do not depend on X, one readily gets 

n t /2 


b(z,z' ;q) = 2zz’ 


d0 *^(8)78 - 2RR* 


(o+a ' ) /2 


tt/ 2 
( a + a *) /2 


d6 cos 20 K x ( S ) / S . 


(4.27) 

In the diagonal case, Eqs. (4.25) and (4.27) simplify respectively to 


2 2 
B(z;q) = 2z d0 K Q (2Rsin 6) - 2R 


ftr/2 


d0 cos 29 Kg(2Rsin 0) , 

(4.28) 


and 


b(z;q) = (z /R) 
fir/2 

- R 


fir/2 


d6 K 1 (2Rsin 0)/sin 0 


d0 cos 29 K^(2Rsin 0)/sin 0 


(4.29) 


If we let here q + 0 we obtain compact integral representations for the un- 
shielded shift functions, viz. 

<1 fit / 2 n 


B(z) = 4z 


b(z) = 2z 


d0 sin 0 KQ(2zsin 0) 


fir/2 


d0 sin 0 K^(2zsin 0) 


(4.30) 


(4.31) 


The connection with the closed form expressions reported earlier"^ is pro- 
vided by Nicholson's formula"^ 


I n (z)K v (z) = (-) n (2/ir) 


ir/2 


d0 cos(n+v)0 K (2zcos 0) 
v-n 


(4.32) 



-15- 


valid when n is an integer , and |Re(v-n)| < 1. Corresponding results are 
not likely to hold for the shielded functions, which are represented by in- 
complete integrals with a variable limit. However, the latter are much eas- 
ier to handle numerically than the principal value integrals in Sec. 3, 

5 . Approximate formulae 

In this section we collect various useful approximate expressions of 

the shielded broadening functions, restricting ourselves to the "diagonal" 

case, which is the most important in practice. 

We begin with asymptotic formulae for A and a, valid when z or/and q 

are large. These are obtained simply by substituting the standard expan- 

18 

sions of Kq and into Eqs. (2.10) and (2.12). The leading terms are 

respectively 

9 9 _9p 

A(z;q) * (tt/2)(1+zVR ) R e * , (5.1) 

a(z;q) 'v (tt/ 2) (l-q 2 /2R 2 ) e 2R . (5.2) 

In particular, if one is interested in the asymptotic behavior for p -* °°, 
then one must let both z and q tend to infinity, while keeping their ratio 
z/q = Zp finite. Thus, when z -*■ »: 

A(z$z/z d ) v (tt/2) (3+1/3) z e 2gz , (5.3) 

a(z;z/z D ) % (tt/26) (0+1/B) e 26z , (5.4) 

where g = (1+1/z^ 2 )^. 

On the other hand, for p -* 0 one has 

a(z;z/z D ) = log(0.68z^/z) , for z^ « 1 , (5.5) 
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whereas ACzjz/z^) -> 1. 

The corresponding formulas for B and b may be derived from the integral 
representations given in Sec. 4. Let us assume first that only z . -*■ <», while 
q remains finite. After changing the variable in Eq. (4.29) to u = 2 R sinQ, 
we expand it as follows : 


2 „2 r2R „„2 2 f 2R 

B = — - — du K n (u) + =- 

R J 2q ° 8R J hq 


du u Kq(u) 


__ 2 , 0 2 r 2R . 

, 5R +3z | , 4 Tr , N , 

+ s— du u L (u) + . . . 

128R J 2q 


(5.6) 


This is readily transformed into an asymptotic expansion by developing the 
coefficients in powers of 1/z and extending all the integrations to infin- 
ity. Eventually 


where 


B(z ; q) 'u ^(qj/z + C 3 (q)/z + . . . , z -> 


x (q) - ^ - q 2 ^Ki 1 (2q) + 2q 2 K ;L (2q) + qK Q (2q) 

= (^ " ^8" + V) K1 l (2q) + (t " t - ) q \ 

*(K) 


C 3 (q) 


(2q) 


qK Q (2q) 


with the new function Ki^ defined by 


18 


r 


Ki 1 (x) = du Kq (u) 
- 1 x 


l-i-1 } 


(5.7) 


(5.8) 


(5.9) 


sari Siro 0 •>- q ■sol t ,b;xsd nsnlio arte nO 

Since Ki^(O) =ir/2 9 in the limit q -y 0 one has C^(0) =tt/4 s C 3 (0) = 9ir/32, 
(?. a&jd we recover the* weli-Jujrqwn |^paiis|.o^^fj : f^z) . 
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Similarly , from Eq. (4.29) one obtains 

3 

b (z;q) 'v c^(q)/z + c^(q)/z + ... , z ® (5.10) 

where 

C x (q) = f | + q 2 ^Ki 1 (2q) - qS^Oq) + qK Q (2q) 

c 3 (q) = (li " Hb" ■ V) Ki i (2q) + (f + 

+ (f "4) qK o (2q) (5 - n) 

In the limit q 0 we have c^(0) = ir/4, c^(0) “ 3 tt/ 32, and Eq. (5.10) re- 
duces to the asymptotic expansion of b(z). 

The above procedure clearly breaks down when q -> •», since all the terms 

in Eq. (5.6) are then of the same order. If z is kept finite an asymptotic 
estimate of B is 9 however,, readily obtained from Eq. (4.28) by noticing that 

the integration interval shrinks as q is increased (tt/ 2 - a 'v z/q) * This 
allows us to write 

B(z;q) i (2z 2 -2R 2 cos 2a) K Q (2 R sina) (ir/2-a) 

-v 2zqK 0 (2q) , q » . (5.12) 

Similarly, from Eq. (4.30) we get 

b(z;q) * zK 1 (2q) , q -► » . (5.13) 

Hence, when z << q both B and b become exponentially small and depend lin- 


early on z. 

Let us consider now the case when z and q tend together to infinity so 
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that their ratio z/q = z D remains finite. In this case no shrinkage occurs, 

but the integrals in Eqs. (4.28) and (4.29) may be evaluated asymptotically 

19 

on integrating by parts. Neglecting higher order contributions one finds 
B(z;z/z d ) ^ (z/Zp 2 ) Ki 1 (2z/z r) ), z -* « (5.14) 

and 

b (z ; z/ Zp) 'v (1/2z d ) Kq^z/z^), z -»■ ~ . (5.15) 

At the opposite limit, as p -*• 0, from Eq. (4.28) one gets B(z;z/z^) -> 0 
and Eq. (4.29) yields 

o 

b(z;z/z D ) -> arctan z D - Jsz D /(l+z D ) , z 0 (5.16) 

The last limit is ^ tt / 2 for z^ >> 1, and = z^/2 ^ or z jy K< ^ • 

6. Numerical results 

Numerical results are not presented here in tabular form for A(z,q) 
and a(z,Zp) since the Bessel functions on Eqs. (2.10) and (2.12) are so 
simple to calculate (see Ref. 18, formulae 9.8.5, 9.8.6, 9.8.7 and 9.8.8). 
Tables 1 and 2 show calculated values of B(z,q) and b(z,z^). Where asymp- 
totic forms could not be used, these functions were calculated both from 
the formulae of Eqs. (4.28) and (4.29) and from the Hilbert transform Eq. 
(3.5) [with direct integration of Eq. (4.26) to give b(z,z D )]. The first 
procedure was by far the simpler, however, overall agreement between the 
two methods of better than 2% was obtained. For values of z greater than 

those in the table, sufficient accuracy can be obtained by using the first 

2 

terms of Eqs. (5.7) and (5.10). Notice that B(z,q=0) = ttz [Kq(z)Iq(z)- 
K^(z)I^(z)] and b(z,z D =«) = ir/2 - ttzKq(z)i^(z) are the straight line 
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results. B(z,q>10) and b (z , z^<. 002) are for all effective purposes neg- 
ligible. 

A(z ,q) , B(z,q) and b(z,q) are plotted in Figs. 2, 3 and 4 for var- 
ious values of q(= p/p^ = z/z D ); a(z,q) was not plotted, but it diverges 
logarithmically at small z and q [Eq. (5.5)]. Notice in particular that 
the functions get small rapidly when q > 1. This is expected since Eqs. 
(5.1) (5.2) (5.12) and (5.13) all predict an exponential fall off when 

q >> 1. This rapid cutoff when p p^, certainly to some extent justi- 

2 12 

fies the usual procedure ’ for treating shielding by a cutoff. To 

further test the validity of these cutoff procedures, in Fig. 5 a func- 

12 

tion F(z,q) is plotted. R(z,q) is essentially [b(z) - b(z )] where 

max 

z = p u)/v. Two cases are considered, firstlv p = p_ and secondly 
max max ’ - max 

p = 0.68 p_. In both cases the agreement between F(z,q) and b(z,q) 
max D 

of Fig. 4 is surprisingly good. (Notice, however, that F(z,q) = 0 for 
p > p .) In fact, when p/p_ <0.1 the agreement for the p =0.68 p_ 
case is better than 5%. This value (0.68 p Q ) was chosen as the cutoff 
since its use in [a(z) - a ( z max )] exactly reproduces a(z,Zp) for small 
z and z D [see Eq. (5.5)], as has been noted in Ref. 7. Actually, for 
q = p/Pp < 0.1 the agreement for all values of z between [a(z) - a ( z mqv ) 3 
and a(z,z D ) is again always better than about 4%; and for q < 0.1 for 
both B(z,q) and A(z,q) the difference between these functions and the un- 
shielded ones (q = 0) is also very small (see Figs. 2 and 3). We there- 
fore conclude that the usual cutoff procedures should certainly be ade- 
quate when simplicity is desired and when p/p^ < 0.1 (which is true in 
most cases of physical interest). However, we believe that the functions 
presented here are of even greater utility, especially if p/p^ should 
get large. 
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APPENDIX 

FORTRAN Programs for B(z,q) and b(z,Zp ) 

(Written by U. Palmer, Joint Institute for Laboratory Astrophysics) 

b(z,z D ) 

b(z,z^) is calculated by FUNCTION FOFZZD(Z,ZD) 

A sample FORTRAIN CALL: B = FOFZZD(Z,ZD) 

FUNCTION FOFZZD requires the following seven routines and a data deck 
marked b(z,Zp) 

1. FUNCTION B2(Z) 

2. FUNCTION FACT(N) 

3 . SUBROUTINE KI1 (X , CAY) 

4. SUBROUTINE INTERP(NX,NY,X,Y,Z,R,S ,T) 

5 . SUBROUTINE LAGRANG (NPTS ,X ,Y ,NP ,LAP , ISTART , IEND ,LEND , R , V ) 

6. SUBROUTINE BESSK(X,CKE,EI) 

7. SUBROUTINE BESSI(X,EI) 

8. data deck b(z,Zp) 

B(z,q) 

B.(z,q) is calculated by FUNCTION FOFZQ(ZZ.QQ) 

A sample FORTRAN CALL: B = FOFZQ(Z,Q) 

The following routines are needed by FUNCTION FOFZQ(ZZ,QQ) 

A. SUBROUTINE INTPB 

B. FUNCTION BQASYM(Z,Q) 

C. ROUTINES numbered 2-7 on the list for b(z,Zp) 

D. data deck marked B(Z,Q) 
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Use of FUNCTIONS FOFZQ(ZZ»QQ) and FOFZZD(Z,ZD) 
b(z s z D ) and B(z s q) are calculated by interpolating in tables B(Z,ZD) 
and B(Z,Q) respectively. Accuracy within the B(z,q) table exceeds 2%. 
When b(z,z D ) exceeds 0.01 the accuracy is also better than 2%, between 
0.01 and 0.001 the first non-zero digit is definitely significant and, 
when less than 0.001, b(z,z D > = 0 will be returned. The tables are pro- 
vided by the data decks marked b(z,z D ) and B(z,q). 

The first call to either routine will cause the reading of its own 
data deck. When inserting the routines into an already existing program, 
care must be taken in arranging the data. Therefore, it would be best 
to make dummy calls to FOFZZD and FOFZQ in the beginning of the main pro- 
gram and at the same time placing the data decks at the start of the pro- 
gram's whole data set. 

For example: 

PROGRAM: PROGRAM MAIN 

DUMMY = F0FZQ(2. 0,2.0) causes reading of table B(Z,Q) and 

stores B(2. 0,2.0) into DUMMY 

DUMMI = F0FZZD(2. 0,2.0) causes reading of table b(z,z^) and 

stores b (2. 0,2.0) into DUMMI 

rest of program 

END 

DATA: Table B(Z,Q) 

Table b ( z s z D ) 

rest of data for Program Main 
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PRQGRAM TESTBS ( I k)OUT * OUTPUT > 00l 

DIMENSION ZZZ<25> *000(25) 002 

DUMMY = FOFZO(2.0»2.0) 00? 

DUMMI . F0FZZD( 2 ®0 » 2 * 0 ) 004 

POINT 800 0O5 

1 ° CONTINUE 00 a 

READ 500*nZZZ*NQOO 007 

TF(NZZZ,lT .0) CALL EXIT 008 

READ ?Ol»(ZZZ( I )» I-I »nZZZ) UOo 

READ 501* (Q00 ( I > » I=1»N00Q) 0-)0 

DO 400 J = 1 «N000 Ol 1 

PRINT 8 70 0i2 

870 FORM *T ( IhO ) Of ? 

0 * OOO(J) Oi4 

DO 4 OO 1 — 2 »nZZZ Oi 8 

Z = ZZZ(I) 028 

BZo = F 0 FZq(Z» 0 ) 0i7 

BZZO = F0FZZD(Z»0) Ois 

PRINT 60l »Z»0»BZ0»BZZD Ol R 

400 continue 0?0 

fOO FORMAT ( 2 O 14 ) 021 

8 O 1 F^ORM AT(5E15»0) Or? 

800 FORM aT ( 2 HJ *1 8X » *Z* * 2 3 X »*0 OR ZD*»13X»*B ( Z *0 )*» 2 2 X »*B <Z*ZO)*»//) Oj? 

801 FORMAT ( 4F20.3 ) 0?4 

END C2S 
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C 

c 

c 


c 

c 

c 


1 


FUNCTION FQFZZDUtZD) ^ , 

COMMON/SET BlZO/Bi 38*21 > »*E<3S> *ZEDE< 21 > *ZL !M( 14) »NZ « NZD *»ZLI M 


DIMENSION CX<3)*CEI<3) 

datailswitch a 1 ) 

GO To <lf2)HSWlTrH 
tSWTTrw « 2 


READ 8< Z*ZD> DATA DEC* 


READ 500»NZ*NZr)*NZLIM 
READ 501* (ZEt I ) *1*1 *NZ) 

DO 11° J = 1 »NZD 
READ 5°1 *^EDE ( Jl 
READ 50 1 » < B 1 I * J) * I =1 *A)Z) 

1 1 0 CONTINUE 

print 6 00.(ZEDEUU) *iu = 1*10) 

PRINT a0 1 ,(ZE(IU)»(B(IU.LU)*LU»i,10),IU = 1 »NZ ) 

PRINT 4 O 3 » (ZEDE (IU ) * IU = 11 *NZD> 

print ft 04 *(ZE(IU)»(BUU»LU)*LU = n»NZD>*IU = 1 *NZ) 
READ 5 O 1 , (ZLIM( I > * 1= 1»NZLIM> 

DO 120 I = 1 *5 

Yd) * AtOGlOIZEDEl 1+16 } > 

12° CONTINUE 

CALCULATE B(Z»ZD) 


2 

•5 

10 


?0 


IF(ZD.GE.O.OOl) GO TO lO 
FORZZD * 0.0 


RETURN 

I F ( (ZD.LT .2*01 .AND. (Z.GT.lO.O) ) GO TO 5 

DO lO j = 2 »NZlI M 

J F { ( ZD.LE.ZEDE( J) ).AND.(Z.GT .ZlTM( J) )>5»2° 
CONTINUE 

I F ( ( ZD.LT..1 ) .OR. (Z.GT.10.0) » GO TO 3O 


IFlZD.GT. 20000.0) 00 TO 40 

I F ( (ZD.GE.lO.O) .AND. (ZD.LE.2000.0) ) GOTO 


*0 


C 

C LAG. INTERP 


CALL I NTEROl nZD*NZ*ZEDE»ZF *B»ZD*Z »TFVP * 4 ) 
FOFZZD = TFMP 
GO To 200 


c p<z*zni asym 

c 

30 ZZD = Z/ZD 
TZZd = 2»0*ZZd 
CALL BESS< (TZZDtOoCEI) 

F<0 - CY( 1 ) 

FKl = C<<2> 

IF(ZD.LT.O.l) GO TO 60 
CALL Kll<TZZD»FFll ) 

TZZdSo * TZZd * ZZD 

FOFZZD a ( i./(i.*Z))»IF)Cll*<l.+TZZDSOUTZZnSO*FKl+TZZO*FlCO) 
GO To 200 

C 

C 

4 O FOFZZD * B 2 -IZ) 

60 TO 200 


360 

361 

362 
3 AO 

364. 

365 

366 

367 
36« 
36° 

370 

371 

3~*2 

3^3 

374 

3-75 

3*7A 

377 

370 

37O 

3p0 

3A! 

3*? 

3 P* 

3 <*4 
3** 
3*>A 
3«7 
3P A 

3po 

300 

3>oi 

303 

303 

306 

30^ 

3«A 

3^7 

3Qfl 

300 

6.00 

401 

402 

603 

604 
603 
40A 
40-7 
40? 
609 

410 

411 
41 7 
41 3 
61 6 
61 A 
41 A 
4? 7 
61 P 
410 
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LOG INTERP 

50 CALL IMTERP(«>»NZ»Y*ZE»B(1»17) *ALOG10(ZD)*Z*TEWP»2) 

FoFZZD a TEMP 
GO To 200 

f>0 FOFZZo = Z * F<1 
?00 TF(FOFZZD.LT*O.OOi ) FOFZZD = 0.0 
RETURN 

500 FORMAT { ?0l4) 

501 form AT ( 1 OFF .4 ) 

ftOO FORMAT! 1H1 »*R(Z»ZD)*»//»* ZO**3 X»i0Fi0.3,//,* Z*»//) 

ftOi FORMAT ( F8»4 » lOFiO.«; ) 

4O0 FORMAT ( im »*8(Z*ZD )*♦//»* ZD* ♦ 3X » 1 1F1O. 3 » // , * Z*»//) 

60a Form AT(F8«4»11Fj0.5) 

END 


4o0 

421 

42? 

42-s 

404 

40F 

Aoa 

4?7 

42 *» 
42° 
4-aO 
4o 1 

400 
4 a O 
4’4 
40c 


FUMCT TOM E?<Z) 

DIMENSION C( 3 ) »FI ( 3 > 

PI = 3.141 59Z654 

CALL BES-Sx CZ»c*EI) 

5? = 0.5 * PT - P T * Z * C<i) * p ! ( "M 

RETURN 

END 


43A 

437 

43P 

43 

44 
44 1 
44 0 


FUNCTION FACT ( M > V7* 

DOUBLE F too 

F = 1 .0 loo 

IF(N .GF. 0) GO To 1.0 i->o 

PRINT aOO iaO 

ftOO FORMAT( 1H0,#A|FGATIVF FACTORIAt*) pi 

call EXIT 140 

lO FACT a 1.0 ip, 

IF(N.GT. I ) . GOTO 2° 134 

RETURN in a 

?0 DO 3° I = T*N tba 

F « F * T 140 

2O CONTINUE lap 

FACT = F 740 

RETURN 100 

END 107 


o o 
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SUBPOUT INE (CH(X.CAY) 1?1 

DIMENSION Elf (7) 7.?? 

DATAtpI =3 *141592 654) ♦( Ek=1 . 25331 414 , 0. 11190289 *0.02576646* 133 

1 0.00933994 *0.00417454, 0.00 163271, 0.00033034) 134 

DATA (EULER* 0.5772156649) 13* 

PI 2 = PI/2.0 136 

if(x.ne.o.o) go to io 1.37 

C ay * PT2 133 

RETURN 1 3° 

10 TF( X.GT .7.0) GO To 200 140 

XT = X/2.0 141 

EPS = 1.0E-9 142 

COEF = -(EUtER+ALOGIXT) )*X 143 

SUM = 0.0 S SUMAsO.O 144 

SUMB * 0.0 5 SUMC = 0.0 145 

|f B = 0 146 

00 1.00 K = KB » 1 00 147 

Tm s 2*K+1 148 

Tie = FACT ( K)**2 149 

T n = XT#* ( 2*Y ) 1 *0 

TA = Tn/(T<*Tm> 151 

TB = TN / (Tk*Tm**?) 152 

SUM = SUM + 1.0 / (<+1.0 ) 153 

TC = (XT**<2*(K>1) )/(FACT(K+l)**2*t2* ( K+l ) + l)) >*SUM 154 

SUmA = SUM a + Ta 15* 

SUMB = SUMB + Tg -»56 

SUMC = SUMC + T C 157 

I F ( ABS (TA) .GT. EPS *ApS (SUmA ) ) GO To 1 00 168 

IFIABS ( TB ) .GT. EPS *ABS ( SUMB ) ) GO To 1 00 15° 

I F ( ABS <TC> »GT . EPS *ABS (SUMf)) GO TO lOO I 4 O 

GO TO 150 161 

lOO CONT INUE 16? 

0 PRTN t 6°0 *Ta *Tb »Tc »SUm»SUma *SUMB *SUmc»K 16? 

fiOO pormaT (40X*7E13.5»I5) 164 

C&LL EXIT 165 

140 CAY a PI? - COEF*SUma - X*SUVB -X *SUVC 16* 

c PRINT 600 »TA»Tb*TC *SUM»SUMA*SUMP*SUmc»K 167 

RETURN 168 

200 X 7 = X/ 7,0 $ SUMO = 0.0 16° 

DO 250 M=1 *7 170 

SUMD = SUMP + (-1 !.0)**(M-i ) * E<(M) / X7**(M-j) PI 

250 CONTINUE 172 

CAY = SUMD/ ( SORT ( X ) * EXP(X)) 17? 

RETURN 174 

END 175 
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SURROUTTNE INTERP(NX»NY»X»Y»Z,R,S»T#NP5 2^5 

DIMENSION X(NX) »Y?NY) sZ(NY»NX) 27A 

DIMENSION DIN« 20) *STO»E< 20) 277 

NPTS = 4 2^® 

r>n i o n = 7*nX 27P 

IF(R.GT,X(N) ) 10.20 2«0 

10 CONTINUE 2«1 

15 PRINT 600.R.S.XINX) »Y(NY) 2*2 

CALL EXIT 2*? 

6,00 FORMAT(*0R»S»X(NX) »Y(NY)*,4E20.9»//) ?«A 

?0 Np = n-i 7<?* 

DO 3° N * 2»NY 

IFtS.GT.Y(N) ) 50*40 787 

3° CONTINUE 7«e 

GO TO 15 7«o 

40 NS = N - 1 2oO 

IF<nS,EQ.1) ns = NS + 1 ?91 

IF(NR.EQ.1> NR = NR + 1 797 

IFImS.EQ. (NY-1 ) ) NS = NS - 1 ?97 

!F(NR.EQ.< NX-1) ) NR = NR - 1 794 

[>n lOO I = 1*4 705 

II = NR-2+I 79* 

DO 5 O J = 1*4 7p7 

JJ = NS-2+J 298 

D I N ( J ) * Z(jj»TT) 790 

50 CONTINUE 300 

CALL LAGRANGU.Y ( NS-l ) »DI N »NPTS .NPTS-l » 1 . 1 » LEND * S »STORE<I) ) 30?. 

100 CONTINUE 302 

CALL LAGRANG<4»X ( NR-1 )» STORE »NoTS » NPTS-1 » 1 * 1 »LFND»R * T ) 3O7 

RETURN 304 

END 305 
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SUBROUTINE LAGRANG (NPTS»X»Y»NP*LAP»I START » TEND # LEND »R »V ) ?Oft 

DIMENSION X(NPTS) *Y(NPTS) *R(IEND) »V( IENDI *DN( 15) *DO< l«i) 307 

DO lOO i * 1*15 30A 

DnU ) S !• 30Q 

1 00 DO ( I ) = 1* 31 0 

LT = NP/? - 1 31 1 

NDMLAP = NP - LAP 31? 

HI s 1 317 

NE » NP 314 

NTemp = 0 31 5 

IR = ISTART 31ft 

1O7 I F ( R ( IR ) - X(NE-LT) ) l02*.i02»l0l 317 

101 IF( (NE. EQ. NPTS). AND. (R( IR).LE.X (NE) ) >G0 TO 102 31 9 

NT ■ NI + NPMLAP 31® 

NE = ME + NPMLAP 32° 

mtemp = ni - i 321 

I F ( ME - NPTS )103»103*104 32? 

104 LEND = IP - 1 323 

PE t Upm 324 

1 0? DO 1 lO x = 1 »N® 3?ft 

xx = x + ntemp ??ft 

DO 110 I = 1 »NP 377 

II = I + NTEMP 37P 

IF(K-T ) I03*ll0*i0ft 390 

lOp DD ( X ) = DO ( X ) « ( X ( XX ) - X(II) > 33O 

11 0 comTImUF 331 

11? V(IP) = 0,0 332 

Do 1 ?0 IT a Nl»NE 337 

IF(R ( IR • - X( IT) >l20*m*120 i»A 

111 V f JR } a Y( I T > 77*; 

GO TO 140 77ft 

1?0 roNT TMUE 777 

DO 130 X * 1 *NP lid 

XX * X + NTEMP 770 

DO 140 I = i»MP 740 

11=1+ NTEMP 741 

IF(X-I)141 *140*141 742 

141 DN(X) = DM ( X ) *(R(IR) - X(II)) 747 

140 continue 7*4 

V(IRi = V ( IR 1 + ( DM*X) * Y ( XX ) X DO ( X ) ) 7*e 

130 CONTINUE 74ft 

149 IF( IR - TEND) 150* 151*151 74 7 

150 JR = IP + 1 74 p 

DO 1?0 mZ a l *NP 740 

170 DN(MZ) « 1,0 7 c 0 

IF( R ( I R ) - X ( NE - LT) ) 1 12* 11?* 161 7 a 1 

161 I F ( NE .LT. NPTS) GO To 162 3ft? 

IF ( ® ( T R ) .LE. X(MPTS) ) GO To 112 7ft7 

16? DO 180 MZ = 1 »NP 7ft 4 

180 DD<MZ) a J 9 0 7ft*; 

GO TO 1 0 1 7 ft ft 

151 LEND * TEND ?ft7 

RETURN 3ft p 

END 3 «P 



- 36 - 

SUB ROUT IMF BESS* IX»CKF*EI) 102 

DIMENSION FIRSK4) »EH 3 ) *C0EF(4) »CKE( *A(iO»4> 

DATA (A = 0. 42278420 * ,23069786* .03488690* 194 

1 .00-563604, , 00010-740, .OOOOOtaO , 3(0,0)* 6.0* 1 on 

2 » 1 8443 144 » -.672788-79, -.18186897* loft 

3 -.0191940?* -.00110404* -.00004686* 3(0.0)* 6.0* 1®7 

4 -.07832388* .02189668* -.01062446* 108 

8 .00587872* -.00251640* .00053208* 3(P.O)* 6.0, 1«9 

6 .23498619* -.03655620, .01504268* ?O0 

7 -.00780353* .003256! a.* -.000«,e?A5* 3(0.0), 6*0 ) 201 

TALL opSS I ( X » F I ) 70? 

T F ( X .Li. ? .0 ) 10,70 777 

1 0 T = X / 7.0 704 

Xp = ALOG(T) 70c 

FIpST(I) = -XP * EI(1) - 0.6?7?i46 6 ?06 

PTDSK?) = X * XP * EI(7) + 1,0 707 

FACTOR a T * T ?Op 

COFF(l) = 1.0 700 

rnEF(?) = 1.0 / x ?*o 

JJ =1 211 

On To 50 7i7 

?0 T= 2.0 / X 7i7 

FIRST ( 3 ) a FIRST ( 4 ) =1.25331414 214 

JJ = 3 316 

C0FF(3) = C 0 EF( 4 ) a 1.0 / (SORT (X) * FXP (X) > ?i 6 

FACTOR = T ? 1 7 

60 JFMD * JJ + 1 21 « 

T = 0 7io 

DO 200 j = J J » JENp 27 0 

I = T + 1 271 

PROD = 1.0 772 

SUM a 0.0 7?7 

KENO = A(l0,J) + O.OOOOOi 774 

DO 1°0 K = 1»KEND 274 

PROP = PPOO * FACTOR 276 

SUM a SUM + PROD * a(K*J) 777 

100 CONTINUE 770 

C<EtI) a COEFlj) * (FTPST(J) + SUM) 770 

200 CONTINUE 230 

0KE<3) = (2.0/X) * CK E ( 2 ) + CKE('l) 771 

RFTUPN 232 

END 237 
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SUBROUTTNE RFSSl (X*EI ) 7*4 

DIMENSION A(l0»4) »FTRST(4) »C0EF(4) »EK *) 7 *« 

OATA(EIpST = 1 . 0, 0. 5 , ?(0, 3989422**) ) * 3*6 

1 (A= 3.5156229* 3.0R99474, 1.206749?. 2*7 

2 .2659732* .0360768* .0045813* 3(0.0), 6.0, 2*8 

3 .87890594, .51498869, .15084934, 2*9 

4 .02658733* . 00301532 * .00032411* 3(0,0), 6.0, 240 

5 .01328592* .00725319,-. 00157565* 241 

6 .00916281* -.02057706, .02635537, 24 2 

7 -.01647633* .00392377, 0.0 , g.O, 24 * 

8 -.03988024, -.00362018, .00163801, 244 

9 -.01031555* ,02282967*-. 02895312* 245 

I .01787654* -.00420059, 0.0 , 8.0 ) ?46 

347 

T = X / 3.75 248 

COEF(l) * 1.0 740 

C0EF(2) = X 7«0 

COEF ( 3 ) = COEF { 4 ) r EXP (X) / SopT (X) 781 

I F ( X .LT. 3.75) 10*20 7*7 

1° FACTOR a T * T 7C7 

JJ = 1 284 

GO TO 50 2*6 

?0 FACTOR * 1.0 / T 7C8 

JJ = 3 247 

50 JEND = JJ + 1 74 « 

1=0 ?eo 

D 0 2°0 J = JJ » JEND 260 

1 = 1+1 261 

PROD = 1.0 767 

SUM = 0,0 76 * 

KEND = A<10»J) + .000001 764 

DO 1°° K = l » KEND 744 

PROD = PROD * FACTOR 766 

SUM = SUM + PROD * A(K*J) 267 

1.00 CONTINUE 768 

EI(I) = COFF(J) * ( FIRST ( J) + SUM ) 340 

700 CONTINUE 7*0 

771 

El ( 3 ) = (-2.0/X 1 * E!( 7 ) + El ( l ) 777 

RETURN 77* 

END 274 
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DATA DECK FOR b (z , z D ) 


3S 21 
* 0 

* 

O >-« 
O > 
-s 

. 00 ? . 

003 

,004 * 

00 * . 

006 

.007 . 

00 @ 

.009 

01 

02 

*03 • 

04 

.05 

06 

07 

.08 

09 

.1 

? 

3 

#4 . 

5 

,6 , 

7 

8 

,9 1 

# 

7. 

# 4 

» 

5. 6 


7 . 8 

• 9 

© 

10 , 



«OOtO 


*0001 

0.0000 

0.0000 

0 » 0000 

0,0000 

0 © npnp 

0 • 0000 

A*pppn 

0«00ap 

o.ooon 

0.0000 

0 .noon 

O.nooo 

o.onoo 

0 , 0000 

p • O^p 

0 ® noon 

P # a A p n 

a # nPnp 

0.0000 

0.0000 

0.0000 

o.nooo 

o.onnn 

0,0000 

n ©noon 

n* pppn 

a* PPO n 

0 #nOOC 

0 «O 0 O 0 

*0050 

*0025 

0.0000 

.0022 

0.0000 

.0017 

0.0000 

.0013 

o.onoo 

.0010 

o.onoo 

.0007 

0 © OOnp 

©0005 

0 # poor 
©0004 

.030 ^ 

©no oz 

*0001 

0.0000 

0.0000 

0.0000 

O.OOO'' 

nonu 

n *oonn 

o*0 onn 

9.0000 

0 9 00(50 

0*0000 

0.0000 

0.0000 

0,0000 

O.OOon 

0,0000 

0 • no am 

U* oonn 

n.OOOo 

OftOOOC 

0*0000 

*01 

ftoo^o 

o.nooo 
. 0 O 4 8 

0.0000 

.0044 

0.0000 

.0039 

O.nooo 

.UU34 

o.onoo 

.0030 

n * on am 
*002 6 

n * noon 
*00 22 

.0019 

«00i 6 

® 00 | 4 

*0092 

0.0000 

0.0000 

0.0000 

0. uouu 

u « OOnu 

V* U(HJO 

U.OOOo 

0 • 0 Or 0 

0®nonn 

o.nooo 

0.0000 

o,nnoo 

0 . 0000 

n, nn^o 

n » nnoA 

0 « n 0 n a 

o. none 

PftpnQC 

0*0000 

0.0000 

0.0000 

0,0000 

O.Oono 

0.0000 

o©oo no 

<>•0000 

.02 

* 0 x 00 

.0099 

.0096 

* 00 92 


*0083 

• u u 7 8 

• UU 7 4 

• UU 69 

• 006* 

© 0060 

*0028 

.0012 

*0005 

©0002 

• UOOl 

u # uonu 

U*UUUU 

U* t>UUfl 

O. 0 O 0 O 

Oe COOP 

o.onoo 

o.nooo 

0 *nonn 

o© 0000 

n« on no 

0 • noon 

n « non?? 

I.UOUUC 

O.OOoo 

0*0000 

0.0000 

0.0000 

0*0000 

0*0000 

U.0000 

0 • uum» 



* 0*500 

.0250 

,0250 

.0248 

.0246 

.0243 

®U240 

*0236 

*0232 

.0228 

• 0224 

*0219 

.0173 

.0131 

.0097 

.0071 

*0051 

©0036 

® OU 26 

.0018 

.0013 

o.onon 

0.0000 

0,0000 

0,0000 

0.0OO r - 

a # n 0 0 0 

n © oonn 

0*0000 

O.OUUQ 

0.9UQO 

0.0000 

.1 

.0502 

0.0000 

.0505 

0,0000 

.0504 

0.0000 

,0503 

0, 00^0 
.0501 

P • on no 
©0499 

C *0C. no 
©U4Q6 

0*0000 
• U493 

.0490 

.0487 

*0483 

.0442 

.0396 

,0350 

.0306 

*0265 

©U229 

* 0 197 

. U 1 6 9 

.0144 

«o n 

e 0 0 0 

*0001 

0.0009 

O.IKK’U 

0*0000 

u • uunn 

U* UUU!.» 

1 1 » u U U 1 ' 

0.0000 

0*0000 

0.0000 

0.0000 

0 • 0000 

O.UOMU 

y u \ *uou 

U iUUIMf 

u « uuun 


.1500 

*0755 

.0764 

*0763 

,0762 

.076^ 

.0759 

.0757 

» f J 7 5 4 

• U 752 

.0749 

*0746 

.0711 

.0668 

,0622 

, 0574 

,0528 

,04S3 

#0 44 1 

.0401 

.0364 

.0121 

.0037 

.0011 

.0009 

.0001 

0. 0000 

o.nooo 

o#nono 

0.000'' 

0,0000 

0*0000 
* 2 

«1°12 

0.0000 

«l008 

0.0000 

•1°°7 

0,0000 
• lOUfi 

O.Oono 
. lOO? 

°« COOO 
. 1003 

0 .0000 

•1*302 

0*0000 
« 1UUU 

*0998 

•0995 

.0903 

*0962 

• O 922 

.0878 

• 0831 

.0783 

.0736 

*0689 

»u643 

.0599 

.0276 

.0118 

.0048 

.0019 

.0008 

.0003 

,0001 

u * 0000 

U.UO'JO 

U.UU90 

0*0000 

0.0000 

0,0000 

0,0000 

0.0000 

0.0000 

0 .0000 

u*uuun 



*3000 

.1538 

.1531 

.1531 

, 1530 

.1529 

.1528 

.1526 

.1525 

.1523 

.1521 

.151.9 

.1493 

.1459 

.1419 

.1375 

.1329 

,12 91 

.1232 

* 1 1 83 

.1134 

*0703 

.0410 

• 0230 

,0127 

.0069 

.0037 

,0020 

«00ll 

• UUU£ 

0.0000 

0.0000 

« 5 
.2637 

0.0000 

*2626 

0.0000 

.2626 

0.0000 

,2625 

o.uoou 

,2625 

0* ouou 
,2624 

0 . uuou 

.2623 

'■) » uuuu 

.2621 

ft 2620 

,2619 

.2617 

.2596 

.2567 

. ’532 

.2493 

«24=,0 

. 240 5 

.2357 

#230*7 

,2256 

« 1738 

.1280 

.0918 

,0647 

e U45 1 

.0312 

»u 2 14 

»<J 146 

®uu99 

.0002 

0.0000 

O.OOno 

0.0000 

0,0000 

0.0000 

0.0000 

0 .0009 

0. oonn 



.7500 

*4035 

.4026 

,4025 

.4025 

• 40?4 

.4023 

#4022 

*4021 

*40 20 

• 40l9 

*4018 

.3999 

.3973 

.3941 

,3904 

.3862 

.3818 

.3770 

.3721 

.3669 

*3104 

.2537 

.2031 

.1604 

.1254 

.0974 

.0752 

.0578 

.04^3 

«UU2 8 

*0002 

0.0000 

0,0000 

0,0000 

0,0000 

n, on no 

n.oooo 

n. oonn 
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1 * 


.5354 

.5347 

.5346 

.5346 

.5345 

.5339 

.5321 

.5296 

.5265 

.5228 

.4396 

.3762 

.3164 

.2629 

.2167 

.001] 
2 . g000 

.000] 

0.0000 

0.0000 

0.0000 

.7520 

.7516 

.7516 

.7515 

.7515 

.7509 

.7491 

.7466 

.7435 

.7399 

.6518 

.5804 

.5096 

,4430 

.3826 

.0074 

*0015 

.0003 

.0001 

0.0000 

2 . 
.9072 

.9070 

.9O7O 

.90 70 

.9069 

.9063 

.9046 

.9020 

,8989 

.8952 

. 80 4O 

.7282 

.6513 

.5777 

.5095 

.0] 87 
3.0000 

.0052 

.0015 

.OOO5 

.0002 

1.0990 

1.0992 

1 .0992 

1.0992 

1.0991 

1.0985 

1.0967 

1.0942 

1.0910 

1.0873 

.<59?8 

.9125 

.3299 

.7492 

.6731 

.0463 

.0177 

.0073 

,0032 

.0014 

5 . 

1.2773 

1.2775 

1.2775 

1.2775 

1,2774 

1.2768 

1.2750 

1.2725 

1.2692 

1.2654 

1.1687 

1.0854 

.9988 

.9135 

.8320 

.0947 
. (1 

.0468 

.0251 

.0142 

.0083 

l u . 

1.4216 

1.4218 

1.4218 

1.4217 

1.4217 

1.4211 

1.4193 

1.4167 

1.4135 

1.4U96 

1.3117 

1.2269 

1.1382 

1.0504 

• 9662 

.1615 

1OO, 

.0971 

.0632 

.0433 

. 0 3 06 

1.5558 

1.5559 

1.5558 

1.5558 

1.5557 

1.5551 

1.5533 

1 .5507 

1.5475 

1.5436 

1.4453 

1.3598 

1 .7703 

1.1816 

1.0964 

.2616 

jOOO. 

.1875 

.1454 

.1181 

• U99U 

1.5693 

1.5692 

1.5692 

1.5692 

1.5691 

1.5685 

1.5667 

1.5641 

1.5609 

1.5570 

1.4587 

1.3732 

1.2837 

1.1950 

1.1097 

.2746 

10000. 

.2003 

.1581 

.1307 

.1114 

1.5707 

1.5701 

1.5701 

1.5700 

1.5700 

1.5694 

1.5676 

1.5650 

1.5617 

1.5579 

1 .4595 

1.3741 

1.2845 

1.1958 

1.1106 

.2754 

20000. 

.2012 

.1589 

.1315 

.1123 

1.5707 

1.5701 

1 .5701 

1.5700 

1.570C 

1 .5694 

1.5676 

1.5650 

1.5617 

1.5579 

1.4596 

1.3741 

1.2846 

1.1958 

1.1106 

. 2764 

*2012 

.1589 

.1315 

.1123 

O 

O 

Oi . 

0 2 

0 6 

1 


4 © 6 & 8 $ 


.5345 

.5187 

,1774 

0.0000 

.5344 

.5143 

.1446 

O.OCnn 

,5343 

.5095 

.1174 

9.0000 

.534? 

.5045 

.0950 

,5340 

,4992 

,0105 

.7514 

.7357 

.3287 

0.0000 

.7513 

.7312 

.2813 

0.0000 

.7512 

.7263 

.24UU 

0.0000 

,7511 

.7212 

.2043 

.7510 

.7157 

.0388 

.9068 

.8910 

.4470 

.0001 

.9067 
. 8864 
.3911 
0.0000 

.90 67 

.8814 

.3413 

0.0000 

.9065 

.8761 

.2974 

.9064 

.8705 

,0730 

1 • 0990 
1.0830 
.6026 
,0007 

1 .0900 
1.0783 
. 5381 
.0003 

1.0989 

1.0732 

.4796 

.0001 

1.0988 

1.0678 

.4271 

1.0986 

1.0620 

.1346 

1.2773 

1.2611 

.7558 

.0049 

1.2772 

1.2563 

.6854 

.0030 

1.2771 

1.2512 

.6209 

.0018 

1.2770 

1.2457 

.5622 

1.2769 

1.2398 

.2163 

1.4216 
1.4053 
.8870 
• 0221 

1.4215 

1.4005 

.8134 

.0162 

1.4214 

1.3953 

.7456 

.0121 

1.4213 

1.3897 

.6836 

1.4212 

1.3838 

.3064 

1 .5556 
1.5393 
1 .0160 
.0649 

1.5556 
1.5345 
.94] 1 
.073° 

1.5555 

1.5292 

.8720 

.0652 

1.5554 

1.5236 

.8087 

1.5552 

1.5177 

.4180 

1 .5690 
1.5527 
1.0293 
• 0Q71 

1.5690 

1.5478 

.9544 

.0861 

1.5689 

1.5426 

.8853 

.0772 

1.5688 

1.537u 

.8219 

1.5696 
1.5311 
.431 2 

■ 1.5609 
1.5535 
1.0301 
.0980 

1 .5698 
1.5487 
.9553 
,0869 

1.5697 

1.5435 

.8862 

.0781 

1.5696 

1.5379 

,8228 

1.56®5 

1.5319 

.4320 

1.5699 

1.5535 

I.O3O1 

.0980 

1.5698 

1.5487 

.9553 

.0869 

1.5697 

1.5435 

,8862 

.0781 

1.5696 

1.5379 

.8228 

1.5695 

1,5319 

.432$ 


4 • 6 *8 1® 2® 
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FUMCTlON FOpZq(ZZ»qq> 0?a 

COMMON /DAT/ Z(45) »Q ( 17 ) »B < 45 * 17 ) »NZ*NQ»NPTS 0?7 

COMMON/OR I G/BB( 45 » 17 ) »pl 0?B 

COMMON/NSm/ BNSM< 45*17) O?' 5 

DIMENSION EI(3)*CkE(3> 0?0 

data ilswitch a 1 > 0?1 

DATACPI a 3.141592654) Oa? 

GO TO ( 100,200) .LSW1TCH 0*’ 

100 L SWITCH = 2 °^4 

READ DATA DECK B(Z.Q) AND GENERATE SMOOTHED TABLES FOR B(Z»0) Oa* 

Oa7 

CALL INTPB 0i *B 

Oar 

CALCULATE B(Z»0) 040 

041 

700 TF(QQ.GT.O.O) GO TO 20l 04? 

CALL BESSK fZZ»CKE»EI) 

FOFZQ = PI * ZZ **2 * ( C<E( 1 )*EI ( l) - CKE( 2) * E I C 2 ) ) 0&4 

RETURN O 4 * 

?0i IF(QQ.LT.lO.O) GO TO 10 04 A 

5 FOFZO = 0.0 047 

RETURN °4B 

lO TEMPB a BOASYM(ZZ.OO) O&R 

IFtZZ.LT.Z(i) ) GO TO 5 0*0 

I F ( ZZ, LE.O .0 5 ) GO T030 0*1 

IF(ZZ.LE.25.0) GO To 40 Oft? 

FQFZo = T EM°B O ^ 7 

RETURN °* 4 

*0 CALL INTERPlNQ*NZ.Q*Z.BNSMtOO*ZZ.T E MP.NOTS) Oe* 

TZ a ZZ**? 0* A 

TQ a 2.0 * QQ 0*7 

RTZTO = SORT ( TZ+TO**? ) 0*Q 

CALL BESSk (RTZTQ*CKE*EI ) ° RQ 

DENUM = PI*ZZ*SORT(TZ+(TQ/Pl )** 2 ) * C^EM) OaO 

FOFZq a TEMP * DENUM 061 

RETURN ° 6 ? 

40 CALL INTERP(NQ.NZ ,0»Z.B »OQ»ZZ *TEM° » NPTS ) 0** 

CALL LAGRANG (NZ*Z*BB(1»1)»4»3»1»1»LEND»ZZ*BBB) 0*4 

tb a temp * bbb * u.o*zz/on * tempb o** 

fqfZq = tb o*a 

RETURN Of 7 

END 0*0 
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SUBROUTINE INTPB 0-r« 

COMMON /DAT/ Z ( 45 ) *Q ( 1 ?) *B U5 ♦ ] 7 ) »K'Z »MO» NPTS 0*>o 

COMMON/OR TG/ P<45 *17) *PI 0 Q 0 

COMMON/ NSM/ BNSM(45»17) O 01 

DIMENSION Cf?)»El(3) Op? 

NZ = 45 0a? 

NQ * 17 Op a 

PT = 3*141592654 Op* 

NPTS = 4 Opa 

Op -r 

READ DATA DECK FOR B(Z* 0 ) Opp 

Op o 

r READ 5 O 1 * (Z( I ) » 1=1 »NZ ) OoO 

DO 1° J = l*NO O 01 

READ 50?*Q( j) Oo? 

READ 502* <P( I»J) »T= 1 *NZ) Oop 

1 0 CONT INUE Or 4 

Ooa 

print table b(z»Q) ooa 

0r7 

PRINT 612* *Q< I) * 1=1 »°) Oop 

PRINT 613 Orr 

PRINT 611 » (Z( I ) » ( P( I » J) * j = l»9> » 1 = 1 »NZ ) 100 

PRINT 612*<QU)*I=lO,NO) 10l 

PRINT 613 10? 

DO 1509 I = 1 »NZ 10? 

PRINT 611 » !Z( I ) , (p( I * J) » J = lO*NO) ) 104 

1 5 OR CONTINUE 105 

DO 15 I = 1*NZ 10a 

107 

smooth table b(z»o) fop z *le. o.o* 100 

100 

5( 1*1) a 1.0 1*0 

temp = 4.°*z< n/pi m 

DO 13 I Y = 1 »N0 11 ? 

TBNSM = SORT! Z( I ) **2+( 2.0*0( IY } )**2 ) in 

CALL BESSfc (T5 NSm*C*EI ) ' 11 l 

DENEW = pI*Z( I ) * SqrT(Z( I)** 2 +f 2*°*Q ( IY)/pl 1> v» * 

BmSm«I» 1Y) = p ( I ’ I Y ) /DENEW 1 1 a 

13 CONTINUE 1*7 

lip 

Smooth Table B<Z»Q) POR Z .GT. 0.05 AND z .LE. 25. 0 !io 

1 ?o 

DO 15 J = 2 *NQ 171 

B ( I » J) a lp(I*J)/p(I*l)) /(TEMP*BQASYM(Zm» 0 ( J) )) l** 

15 continue 11 * 

5°1 FORMAT ( lOEg.O ) 174 

5 O 2 FORMAT ( 6E12.4) 1?c 

611 PORMAT { F 1 O .3 »9E14*4 ) 1 -»a 

612 FORMAT ( 1H1 »*B < Z »Q) * » / / *7X »*0* » 2* *9El4 *3* / > 1?? 

613 FORMAT ( m ,* Z** ' ) 1 -> o 

RETURN 1?o 

END i^O 
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FUNCTJON BQASYm(Z»Q) Oao 

DIMENSION CKO)*CEI<3> O"’ 0 

TWOQ = 2.0 * 0 0-J] 

CALL BESS|((TWOO»C|C»CEI) 077 

CALL KI1«TW00»CAY) ° 77 

QSQ = 0**2 

BOASYM s ( j.O/Z )*( (0.5-QSO)*CAY+Q*C<( 1) + 2»°*QSQ * C< ( 2 1 > 07«i 

return 07 * 

END 077 
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DATA DECK FOR B(z,q) 


.001 

002 *003 

.0 04 

.005 .006 

.007 

008 .009 

.Din 

.020 . 

O 3 O *040 

,0 50 

.OgO *070 

»OgO 

090 *100 

« 2 A 0 

®300 » 

4 OO *500 

*600 

.760 .000 

.900 l. 

000 1.700 

1 • <*''0 

1.600 l» 

gOO 2*000 

3.000 4 

.000 5.000 

6*000 8 * 

0 nr ir . 000 

12*000 

14.000 16. 

000 18,000 

20.000 25 

,000 




1 * OOOOf-0 8 
2.0495E-05 

7.3269E-05 

1.5339E-04 

2.5824E-U4 

3. 8597E-04 

5.35199-04 


7.0472E-04 

8 < 9362E-04 

1. 1010E-03 

1.3262E-03 

4. 4350E-03 

8.8763E-03 


1 .4272E-02 

2 *06649-02 

2.7577E-02 

3.5205E-02 

4.3356E-02 

5.1957P-U2 


6.0942E-02 

1*621 4E-0 1 

2.6576E-01 

3 » 5 9 84E-0 1 

4 * 366 2E-0 1 

4.9P0&E-U1 


5.i397E-°l 

5®i594E-0i 

5 • 9594E-03 

6 * 0392 E-0 1 

6 .O 312 E-O 1 

5 . 8 OO 5 E-O 1 


5 .4588E-0 1 

5 • 0684E-0 I 

4.669 uE- 01 

? • 05 A5E-0 1 

2 * 15 8 2E-0 1 

1 *66269-0 1 


1*3578 E- 01 

1 » 0006E-0 1 

7.P469E-02 

6 . 5 9 7g E—0 2 

5. 6479E-02 

4.9306E-02 


4.3785E-02 

3.9379E-0? 

3. 147OE-0? 





6 » 00 OOe-O 3 
4 .8939E-05 

1*076 1 E-0 4 

1 . 8318E-04 

? . 79 5 1 E-0 4 

3. 97919-04 

5.383 2 c -0u 


7.0001E-04 

B » 8208E-04 

3 .0836E-03 

3 .3037E-02 

4.3862E-03 

8 • 7779E-03 


1.4209E-02 

2 * 0498E — 0 ? 

2.751 0E-02 

3.5136E-02 

4. 3286E-02 

5.1887E-02 


6.0871E-02 

1 • 620 7E—0 1 

2 » 6569E-01 

3 • 5 878E-0 1 

4.3656E-01 

4.980lE-0i 


5*4391 E-0 1 

5 . 7590E-01 

5 . 95 89E-0i 

6.058 8E-0 1 

6 *09099-0 1 

5 « 8 OO 2 E-O 1 


5.4586E-01 

5.0682E-01 

4 . 6688E-01 

3.0544E-01 

2 • 1 581E-01 

1.6624E-01 


1 *357RE-0i 

1 » OOO 5 E-O l 

7, 1466E-02 

6 » 5 9 76 E-0 7 

5 « 647‘ 7 E-02 

4 . 03 04 E —0 2 


4.3 784E— 0 2 

3.9379E-02 

3. 1471E-02 





1 . OOOO 9 -O 2 
8.1654E-05 

1 .6695E-04 

2 . 6092E-04 

3.6659E-04 

4. 8655E-04 

6 . 2 26 5 E— Of. 


7.7606E-04 

9.4739E-04 

1 . 1368E-0? 

3 .3444 E—0 3 

4. 3253E-03 

8 • 6646 E-O 3 


1 .4067E-02 

7 * 0338E-0? 

2 * 7338E-07 

7.4966E-07 

4.3101 E—0 2 

5.1 6 c, 7E-0 •? 


6 » 0 678E-0 2 

1* 61879-01 

2.6550E-03 

3*58 5 9 E— 0 1 

4.36 9 9E-*0 1 

4.9785E-01 


5*43 77E-0 1 

5.7676E-03 

5.9577E-01 

6.0576E-03 

6.0298E-01 

5.7=>95E-01 


5*4 5 789-0 1 

5*06' 7 5F-0i 

4 . 66 82^-0 1 

3 . 0 5 tO E- 0 1 

2 . 1 578E-03 

1 • 6622E-0 1 


1 .3576E-01 

1.0004E-01 

7 • 9 4 5 4 F - 0 2 

6.5965E-U2 

5. 641SE-02 

4 . 9296 E—0 7 


4.3776E-02 

3.9372E-02 

3. 1463E-02 





5 .OOOOE-O 2 

2 .4279E-04 

4, 8608E-C4 

7. 3035E-C4 

9.761 PE-U4 

3 . 2240E-03 

3 .4742E-03 


1 .7275E-03 

1.9841E-03 

2 . 2446E-03 

2 . 5u 93E-u 3 

5. 4 76 2 E-0 3 

9.21 IIE-O 3 


1 » 3 840E-0 2 

1.9355E-02 

2.5681E-02 

3 » 2 7 ZOE-u 2 

4.0082 E—0 2 

4.8571E-07 


5 * 7? 1 1 E-O 2 

1 » 5688 E-O 3 

2. 6C24F-01 

3 . 534 3 E-0 1. 

4 . 3144E-01 

4.9318E-01 


5 » 3939E-0] 

5*71 67E-C 1 

5.91 95E-0 1 

6.0221 E-0 l 

5. 99?9E-0l 

5 . 772 3E-0 1 


5 . 434 OE-O l 

5 » 0464E— 0 ] 

4 • 6 4 9 4 E - 0 1 

3 • U4 20’E-U 1 

2. 1490E-01 

1 .6553E-01 


1.3519E-01 

9.9613E-02 

7.9116E-02 

6.5684E-02 

5.6177E-02 

4.9086E-02 


4.3590E-02 

3»9204E-02 

?. 3 33 0E-O2 





1 , OOOOe-Oi 
3 • 5 072E-04 

7*01 31E-04 

1 .O624E-03 

1 .4040E-03 

3 .7563E-03 

2.1094E-03 


2*4636^-03 

2 9 81 90E-0 3 

3 . 1 757 E-O 3 

3.5339E-03 

7 * 2351E-03 

1. 1253E-02 


1 * 57 O 4 E-O? 

2 9 O 668 E-O 2 

2* 6 186- “0 2 

3*227uE-02 

3 . 8906E-02 

4 . 6 O 6 IE-O 2 


5.3697E-02 

1 * 4668E-0 1 

2 • 4764E-0 1 

3. 40 18 E-0 1 

4* 1830E-01 

4.8050E-01 


5.2733F-01 

5.6030E-01 

6. 8 1 2 8E-0 1 

c .9222E-0l 

5.913 6E-01 

5.69595-01 


5.366TE-01 

4.9868E-01 

4. 5Q61E-01 

3 . 00 8 3 E- 0 1 

2. 1245E-01 

1 .6360E-01 


1*33599-03 

9.8420E-02 

7*81 64E-02 

6.4892 E-0 ? 

5 » 5498E-02 

4,8493 E-02 


4.3060E-02 

3.8726E-02 

3.0945E-02 





1 , 5 OOOE-Oi 
4*11 75E-04 

8*2359E-0i 

1 » 2356E-U? 

1 « 64 7 9E-U3 

2.0605E-03 

7 * 47359-03 


2.8870E-03 

3.3012E-0? 

3.7160E-O3 

4*13 16E-U 3 

8.3473F-03 

1,27269-02. 


] .7340E-02 

2.2246E-07 

2.-7490E-U2 

3*31uZE-u2 
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